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IMAGE  UNDERSTANDING  AND  INFORMATION  EXTRACTION 


ABSTRACT 


V 

\ 


This  report  summarizes  the  results  of  our  research  program  on  Image  Under¬ 
standing  and  Information  Extraction  supported  by  the  Defense  Advanced  Research 
Projects  Agency  under  Contract  F30602-75"C-0150.  The  report  covers  the  period 
February  1  to  April  30,  1976. 

-''The  objective  of  our  research  Is  to  achieve  a  better  understanding  of 
image' structure  and  to  use  this  knowledge  to  develop  techniques  for  Image 
analysts  and  processing  tasks,  especially  Information  extraction.  Our  emphasis 
Is  on  syntactic  decomposition  and  recognition  of  imagery  based  on  scene  analysis. 
It  Is  our  hope  that  the  results  of  this  research  will  form  the  basis  for  the 

development  of  technology  relevant  to  military  applications  of  machine  ex- 

/ 

traction  of  Information  from  aircraft  and  satellite  imagery. 


\ 


\ 

\ 
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IMAGE  UNDERSTANDING  AND  INFORMATION  EXTRACTION 
Research  Summary 

This  report  summarizes  our  research  progress  during  the  period  February 
1  to  April  30,  1976,  In  Image  Understanding  and  Information  Extraction. 

Our  research  objective  Is  to  achieve  a  better  understanding  of  image 
structure  and  to  use  this  knowledge  to  develop  techniques  for  Information  ex¬ 
traction  from  Imagery.  It  Is  our  hope  that  the  results  of  this  research  will 
form  the  basis  for  the  development  of  technology  relevant  to  milstary  applica¬ 
tions  of  machine  extraction  of  information  from  aircraft  and  satellite  Imagery. 

Our  research  projects  fall  Into  five  heavily  overlapping  areas:  Image 
Segmentation,  Image  Attributes,  Image  Structure,  Image  Recognition  Techniques, 
Preprocessing,  and  Applications. 

IMAGE  SEGMENTATION  -  We  pursue  two  approaches  to  Image  segmentation: 
edge  detection,  and  region  growing.  In  edge  detection,  the  thrust  of  our 
research  Is  to  use  syntactic  methods  to  help  edge  detection  in  noisy  and  blur¬ 
red  Imagery.  As  a  step  along  that  direction,  we  are  studying  characteristics 
of  digital  edges.  Almost  all  past  works  on  digital  edges  and  curves  dealt 
with  the  Ideal  situation,  while  we  are  interested  mainly  In  real-life  digi¬ 
tized  Images.  Our  experimental  results,  reported  by  Tang  and  Huang,  indicate 
that  the  properties  of  real-life  digital  straight  edges  are  quite  different 
from  those  of  ideal  edges.  In  the  same  report,  a  technique  of  recognizing 
real-line  digital  straight  edges  is  also  presented. 

A  by-product  of  our  edge  detection  research  is  a  technique  of  accurately 
estimating  edge  locations  which  holds  great  promise  in  mensuration  applica¬ 
tions.  This  technique,  reported  by  Burnett  and  Huang,  Is  based  on  the  Viterbl 
algorithm.  It  can  take  account  of  arbitrary  edge  profiles  and  film  noise  and 
Is  computationally  simple. 
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In  region  growing,  we  are  studying  both  region  merging  techniques  and 
different  similarity  measures  among  regions. 

IMAGE  ATTRIBUTES  -  We  are  doing  both  texture  and  shape  analysis.  A  new 
class  of  texture  descriptors,  the  max-mln  descriptors,  has  been  developed. 
These  descriptors  are  very  simple  to  compute  and  perform  extremely  well  In 
various  classification  problems.  We  are  currently  exploiting  the  use  of  max- 
mln  descriptors  in  texture  boundary  detection  and  region  growing. 

In  shape  analysis,  we  have  done  extensive  study  on  Fourier  boundary 
descriptors  and  are  striving  to  Improve  their  performance.  We  are  also  devel¬ 
oping  grammars  for  various  classes  of  objects  of  military  significance,  such 
as  airports,  and  tanks. 

IMAGE  STRUCTURE  -  The  thrust  of  our  research  In  this  area  Is  to  use 
syntactic  methods  to  do  scene  analysis.  Fu  and  Li  report  on  the  use  of  tree 
grammar  In  detecting  highways  and  rivers  In  LANDSAT  Imagery.  Tree  grammar  was 
also  used  In  helping  scene  segmentation  (Fu  and  Keng). 

IMAGE  RECOGNITION  TECHNIQUES  -  We  pursue  twc.\J:opIcs:  the  use  of  branch 
and  bound  techniques  in  solving  recognition  problems,  and  the  use  of  context 
In  statistical  classification.  In  Swain  and  Yu's  report,  we  see  that  the  use 
of  context  Indeed  increases  the  classification  accuracy.  However,  It  is 
computationally  tedious  on  conventional  serial  computer.  We  are  therefore 
looking  Into  the  use  of  the  I II  lac  4  (via  the  ARPANET). 

PREPROCESSING  -  The  aim  of  preprocessing  is  to  change  the  image  to  a 
form  which  Is  more  convenient  for  Information  extraction.  Most  Images  are 
degraded  by  noise  and  blurring.  The  reduction  of  noise  and  the  sharpening  of 
the  Image  generally  facilitate  information  extraction.  Several  noise  reduc¬ 
tion  techniques  are  compared  by  Yoo  and  Huang.  It  was  found  that  the 


synthetic  highs  technique  works  best  in  balancing  the  noise  ieve)  and  the  edge 
sharpness. 

APPLICATIONS  -  We  are  working  on  several  applications,  two  of  which  are 
reported  here.  Wai iace  and  Wjntz  present  some  recent  results  on  using  Fourier 
descriptors  for  airplane  classification.  Hi tchel i  and  Chen  show  results  of 
reducing  cloud  and  haze  in  LANDSAT  imagery  using  three-dimensional  digital 
filters. 


APPLICATION  OF  A  FINITE  STATE  MARKOV 
PROCESS  MODEL  TO  IMAGE  MEASUREMENT 


J.  Burnett  and  T.  S.  Huang 

I •  Introduction 

Our  last  report  [1]  briefly  mentioned  the  possibility  of  using  the 
Vlterbl  algorithm  In  conjunction  with  a  finite  state  Markov  process  model  for 
making  accurate  measurements  from  noisy  and  blurred  Images.  Here  we  expand 
on  this  idea,  derive  the  algorithm  and  formulas  for  calculating  the  perfor¬ 
mance  of  the  algorithm. 

A  Computer  simulation  of  a  specific  photographic  imaging  system  shows 
the  algorithm  to  produce  asymptotically  efficient,  and  unbiased  estimates  of 
object  sizes. 

Maximum  A  Posteriori  Probability  Sequence  Estimation 

The  maximum  a-posterlori  probability  (MAP)  estimate  of  a  sequence 
given  a  sequence  z_  (£  being  a  degraded  and  noisy  version  of  I)  Is  defined  as 

A  A  A  A 

a  sequence  J_  *  (l  j » im)  such  that  P(ljz>)|,/j  Is  a  maximum.  To  calculate 
* 

I  a  model  Is  needed  for  the  relationship  between  J_  and  z.  We  assume  the 
observation  model  shown  in  Fig.  1.  J_  ■  (l j , l2» . .. , 1^)  Is  the  sequence  of 
Ideal  light  Intensities  with  l^  the  light  Intensity  of  the  kth  sample  point 
entering  the  Imaging  system.  Each  i^  can  assume  one  of  G  possible  values 
a|,...,aQ.  For  example  J_  might  represent  the  sequence  of  reflected  light 
Intensities  from  a  scan  line  of  an  aerial  photo  of  a  bridge  across  a  river. 

In  this  case  there  would  be  two  possible  intensity  levels:  a^  corresponding 
to  the  light  reflected  from  the  water  and  a£  corresponding  to  the  light 
reflected  from  concrete  (or  whatever  construction  material  was  used  in  the 
bridge).  The  state  at  position  k,  t^,  is  defined  to  be  a  set  of  adjacent 
Intensities  (I,  . I.,lk  Since  each  I.  can  assume  only  a  finite 
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number  of  values  each  Is  one  of  a  finite  set  [Sj.  .  .  Sp],  Further  (to 
within  boundary  conditions)  there  is  a  one  to  one  correspondence  between  the 
state  sequence  £  and  the  intensity  sequence  J_. 

The  system  h(*)  represents  the  degradation  of  the  sequence  J_.  In  the 
case  of  photographic  imagery  this  includes  blurring  due  to  scattering,  dif¬ 
fraction,  camera  motion,  etc.  as  well  as  the  nonlinear  relationship  between 
light  Intensity  and  film  density.  The  only  assumption  that  we  make  on  h  is 
that  there  is  a  one  to  one  correspondence  between  Y_m  (yj».«.»ym)  (where 

Yk  "  h^nk^  and  H* 

N  Is  a  sequence  of  independent  noise  samples.  We  do  not  rule  out 
dependence  of  the  noise  parameters  on  the  signal,  however.  For  example  film 
grain  noise  is  approximately  normal  with  a  standard  deviation  proportional 
to  the  signal  level. 


The  Algorithm 


By  definition  the  MAP  sequence  estimate  1  of  I  Is 


(1)  P(jJZ)Hf  is  a  maximum 

but  since  there  is  a  one  to  one  correspondence  between  I  and  T^(l)  is 
equivalent  to 


(2)  P (rj_| Z_)  Is  a  maximum 


However, 


(3)  p(nJl>  ■  p(z»a)/p(£)  ■  p(ila)  p(d)/p(z) 


Due  to  the  Markov  assumption  of  and  the  independence  of  the  noise 

M 

p(tl)  -  -T  P(nk+1  In.) 
k-l  K  1  k 
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P(2|a)  -  P(2|Y) 


H  M 

*  p(zjYiJ  "  P (2.  ! h(n. ) ) 

k«l  k  k  k-l  k  K 


Thus  we  want  to  maximize 


H 

(*>  P(*WV 

or  equivalent  minimize 

M  H 

(5)  ^  -  An  P(nA+1|nA)  -  An  p(*£|h(na))  ^y(n£) 


By  assigning  a  cost  or  length  of  r(ri£)  to  each  branch  of  a  trellis  we  can 

A 


see  that  the  MAP  estimate  n_  represents  the  lowest  cost  or  minimum  length  path 

through  the  trellis, 

A 

Let  rij  be  a  sequence  of  states  starting  at  some  initial  state  at  posi¬ 
tion  one  and  ending  at  state  n  at  position  A .  In  general  there  will  be 

£ 

several  possible  paths  (or  sequences  rij)  through  t>e  trellis  that  pass 
through  state  n  at  position  1.  Denote  this  set  of  paths  by  E.  Let  r> |  e  E 

be  one  of  these  sequences  but  with  the  additional  restriction  that 

*A  A  ^  ,A  *  a  A 

r (nt )  ■  II  r(n.)  <  r(n.)  for  any  other  sequence  n,eS.  (If  the  two  paths 

1  j-1  J  1  1 

have  equally  low  cost  any  reasonable  procedure  for  deciding  between  them 

will  do.) 


Thus  Ti|  (called  the  survivor  sequence  or  survivor)  Is  the  minimum  cost 

sequence  or  path  from  the  fixed  Initial  state  to  the  state  r|  at  position  A. 

Now  If  the  minimum  cost  complete  path  from  position  1  to  position  m  passes 

through  state  q  at  position  A  It  must  have  rij  as  Its  Initial  segment  (If  It 

AA 

did  not  we  would  replace  the  Initial  segment  with  and  get  an  even  lower 
cost  path,  a  contradiction). 


Of  course  we  do  not  know  that  the  minimum  cost  complete  path  passes 
through  state  n  at  position  l.  However  we  do  know  that  It  must  pass  through 
one  of  p  possible  states  at  position  £.  Thus  at  any  position  £  we  only  need 

A  q  a  0 

store  at  most  p  survivor  sequences  and  their  costs  r(rij).  To  get  to 
position  l+|  we  need  only  extend  all  position  £  survivors  by  one  unit,  com¬ 
pute  the  costs  of  the  possible  extensions  from 

(6)  r(n|+l)  -  r(n,>  +  r(ni+l) 

and  for  each  of  the  p  possible  states  at  position  £+1  select  the  lowest 
cost  path  ending  in  that  state  as  the  position  £+1  survivor.  In  summary: 

Storage:  £  (position  Index) 

Hj  (p  such  £  point  survivor  sequences) 

£ 

r(nj)  (costs  of  each  of  the  p  survivor  sequences) 

Initialization:  £-0 
*0 

Hq  ■  n0  for  each  possible  Initial  state 

r(nJ  ■  "An it  where  u  Is  the  a-prlor!  probability 
u  n0  n0 

that  the  Initial  state  is  Hq  (If  known) 

If  the  a  priori  probabilities  are  not  known  then  any  reasonable  Initial 
erst  assignment  (such  as  r(nQ) «  0  for  all  possible  Initial  states  n)  will  do 
Recursion 

For  each  of  the  p  possible  states  at  position  £+1  compl  y 

r|,wV 6  r(r>f> +  r<v,) 

for  each  possible  find  r(n]  )  «  min  r(n£+.  .n^)  store  r(n  ')  and  the 


corresponding  survivor  sequence.  At  position  M  there  will  be  at  most  p 
survivor  sequences,  one  survivor  sequence  terminating  at  each  permissible 

A  A 

position  N  state.  Denote  by  r[  the  lowest  cost  of  these  sequences,  jx  is  the 
HAP  sequence  estimate  of  rj. 

Example.  Suppose  It  Is  known  that  at  sample  point  number  one  the  local 

gray  level  Is  zero,  that  at  position  eight  the  local  gray  level  Is  three  and 

that  someplace  between  these  two  points  a  step  change  between  the  levels  zero 

and  three  occurred.  The  degrading  system  Is  linear  and  shift  invariant  with 

iitpulse  response  h-1 *  h^ *  h^  ■  y .  The  noise  Is  white,  Gaussian  with  variance 
2 

o  .  The  observed  sequence  £  ™  (-.5*  +.25,  +.75,  .5*  2.8,  2.7»  3.3,  3.1). 

The  possible  states  are 


S]  -  (0,0,0) 

S2  -  (0,0,3) 

h(Sj)  -  0 

h(S2)  -  1 


S3  -  (0.3,3) 

-  (3,3,3) 

h($3)  -  2 

h(S4)  -  3 


By  assuming  that  all  permissible  changes  of  state  are  equally  likely 

)  -  (z£-h(n))2. 

The  trellis  for  this  example  is  shown  in  Fig.  2.  The  permissible  paths  are 
shown  by  dashed  lines.  The  numbers  are  the  costs  of  the  survivor  sequences  to 
each  state.  The  HAP  estimate  is  shown  by  a  solid  line.  Its  total  cost  Is 

A 

1.955  and  corresponds  to  I  «  (0, 0,0, 0,3, 3, 3, 3) . 

The  minimum  cost  or  minimum  length  path  through  a  trellis  problem  and 
various  solution  have  been  around  for  some  time  [2].  The  algorithm  presented 
above  (commonly  called  the  Vlterbl  algorithm  or  VA)  was  presented  by  ViterbI 
[9]  as  a  technique  for  decoding  convolutional  codes.  The  algorithm  has  since 


-An  Pfr^j In^)  terms  are  the  same  and  can  be  Ignored.  Thus  T(n^ 
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been  used  by  Forney  [33  as  a  solution  to  the  Intersymbol  interference  problem. 

The  algorithm  also  has  applications  In  text  recognition  [4]  since  it  can 
exploit  Harkov  dependence  In  English  text  and  can  be  applied  to  scene  analysis 
problems  [5]  due  to  Markov  relationships  between  objects  in  scenes. 

The  VA  can  be  used  to  make  measurements  of  objects  In  digitized  images 
providing  the  object  whose  length  or  width  ts  to  be  measured  is  distinguished 
from  Its  background  by  abrupt  changes  in  reflected  light  intensity  at  Its 
edges.  A  HAP  sequence  estimate  of  a  scan  line  across  the  object  can  be  cal** 
culated  and  the  edge  locations  relative  to  the  start  of  estimated  line  can  be 
subtracted  to  produce  an  estimate  of  the  size. 

In  the  next  sections  formulas  will  be  derived  for  predicting  the  per** 
formance  of  the  VA  at  locating  edges  and  estimating  widths. 

Error  Analysis  for  Step  Edges 

In  the  previous  section  we  presented  the  VA  as  a  possible  technique  for  making 
measurements  from  digitized  Images.  The  measurement  technique  consists  of 
making  a  HAP  estimate  of  a  scan  line  across  the  object  of  interest  and  sub¬ 
tracting  the  position  of  the  boundary  points  of  the  estimate  of  the  scan  line 
to  get  an  estimate  of  the  width.  In  this  section  we  examine  the  accuracy  of 
this  technique.  In  particular  we  will  calculate  the  probability  of  mlslocat- 
fng  a  step  edge  by  |n [ (n «  +1, +2,...)  points.  It  is  assumed  that  a  change  In 
brightness  from  level  a^  at  some  initial  position  to  level  a^  at  position  m 
has  been  detected.  The  observed  signal  is  blurred  and  noisy  so  that  the 
location  of  the  step  change  between  levels  a^  and  a^  is  uncertain. 

A  A 

Define  an  error  event  En  by  the  conditions  an<*  *  1  n  | + 1  *"  ^  k+ 1  n  |  -*■  1 

but  lj  j*  ij  for  k+l  <,  J  <.  k+|n|.  If  p  Is  the  number  of  states  then  there  must 

* 

|  n  |  +  p-2  (|n|  >_  !)  state  disagreements  (see  Fig.  3  for  p»A,  n* -2) ..  If  £ 

A 

Is  the  output  vector  corresponding  to  J_  and  ^  Is  the  vector  corresponding 
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to  the  t.  ue  sequence  J_  then  ^  w(1 1  be  decided  over  ^  if 

k+|nl+t  a  k+|n|+t 

(7)  2  An  P(z,|y.)  >  l  An  P(z.|y.) 

j-k-t  J  J  J-k-t  J  J 

where  t  -  f(P'1)/2]  -  largest  integer  (P-1 )/2 

2 

If  the  noise  Is  normally  distributed  with  zero  mean  and  variance  a  then  (7) 


becomes 

k+|n|+t 

a  k+|n|+l 

(8) 

1  (VV  ’ <  z 

j-k-t  J  J  j-k-t 

or 

ur-rii2  < 

nr-rti2 

where  Z'»  (Vt’*" '  V“”Zk+  n  +t) 


(yk-t . yk+t+  j  n  | 5 

(yk-t  *  *  *  *  *y  k+t+  J  n | ^ 


The  vectors  z£  Y?  V'can  be  viewed  as  three  points  in  |n|+2t+l  dimensional 

space  that  define  a  plane  (see  Fig.  5). 

Define  <Z,Y>  »  z"*"y  -  E  z.y.  =  Y^Z.  Now  the  condition  for  deciding 

j  J  J 


Y  becomes 


jr^rn  W-zi-*  < 

Ak 

(the  reason  for  the  normalization  by  ||Y-Y||  will  be  obvious  shortly) 

— [<y;y -  <r,Yi&  <  <y^y£z*  -  <r,Y-z$] 

1 1 Y— Y 1 1 

— —  .  [<£r^z^  -  <z;y*  <  <y;y*z*  -  <r;Y?>] 

| Y-Y | | 
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l|Y-Y|| 


[<vj  r-z'>  -  <fx>  “  <x»  l'>  <  <y^  y -  <z;  y->  -  <y;  y^I 


—Ik —  hir  y1!^  -  <y;  Yit->>  <  <y;  y^z">  -  <y;  v*r>] 

ll-ll  I 


A 


(9)  <r-  y;  r-  r> ,  <r-  s  r-  r> 

i ir  -  ?'i i  nr -i'll 


Thus  an  error  will  occur  In  deciding  between  Y*and  Y* whenever  the  distance  from 
^  Y**~  Y* 

Ifto  Y'along  the  —  —  axis  Is  less  than  the  distance  from  Z^to  Yfalong  the 

iiy:-  r  ii 

...Y-Y 

same  axis.  Since  Z^»  Y*+  N"the  projection  of  7?-  Y*on  - * —  Is  determined 

iir-rii 

by  the  projection  of  the  noise  N'  on  the  axis.  This  projection  Is  a  linear 

2 

combination  of  normal  random  variables  with  zero  mean  and  variance  a  and  hence 
Is  also  a  normal  random  variable  with  zero  mean  and  (due  to  the  normalization 

^  A  A 

by  | |Y-Y| |)  variance  cr  .  Therefore  the  probability  that  Y  will  he  decided  over 

the  correct  path  Y  Is  the  probability  that  a  normal  random  variable  exceeds 

half  the  distance  between  Y**  and  Y**.  With  T  ■  1 1 Y**—  Y^H  then 

—  —  n 

/i ps i_  f  v/ *•  1  _*n/2^  . 


(10)  Prob  (decide  Y"*  over  Y'’)  **  f. 


n  /  2  to 


-  • 

A  A 

Now  Y  -  Y  ■  (I  -I)  *  h  (assuming  a  linear  and  shift  Invariant  degradation) 

I  -  *1  -  (0,0,0,A,A,A,0,0) 
where  A  ■  +(a^  -  3j) 
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Let  Uj  (K)  be  the  unit  step  response  of  h.  Then  |  j Y"-  Y''|j2  ■  f  | Y  -  Y||2 
r  r  *  \2  t2 

-  J  W  -  rn 

r?  -  5  [(i,-1.)  *  h]2 

"  J-l  j  J 

0  k+|n|+t  , 

-A 2  Z  [U.(j-K)  -  U.Q-Inl-K)]2 
J-k-t  1  1 

-  |n|+t  , 

(11)  -  A2  I  [U,  (J)  -  U.(j-|n|)]2 

J-t 

For  any  value  of  n  there  are  two  possible  edge  locations  estimates  that 
are  |n|  points  In  error  corresponding  to  location  estimates  that  are  n  points 
before  and  n  points  after  the  correct  location.  Thus  the  probability  of  an 
error  event  for  step  edges  with  known  levels  Is 

(12)  P(En)  -  2Q(^) 


where  T  Is  given  by  (11). 
n 

Thus  far  we  have  only  considered  the  case  of  constant  noise  variance. 
However,  for  some  tyoes  of  noise  such  as  film  grain  the  variance  of  the  noise 
varies  with  the  signal. 

Consider  again  the  problem  of  trying  to  locate  a  step  edge  between  levels 

2 

a,  and  a,  but  now  assume  that  the  noise  variance  of  the  Ath  sample  a  depends 
1  z  y& 

on  y^.  Equation  (7)  becomes 


k+ 1 n  j  +t 


(13)  V  (/2tt  a*  )  -  (z,-y,)2/2a2  >  E  -  in  /2tF  a 


k+ I n I +t 


J«k-t 


YJ 


j  T 


Yj  J-k-t 


*  (Zj-Yj)2/2^ 


J 


or  In  an  obvious  notation 


(H)  nr-rii  °2  <  iii'-  x'i ict^+c- 

^  i 


where  C*  ■  Z  Jin  . 

yy  yJ 

*\ 

The  exact  probability  of  deciding  Y  over  Y  can  be  calculated  from  (14) 
by  a  procedure  given  in  [6)  though  It  Is  very  difficult  and  will  not  be  con- 

A 

sldered  here.  However,  one  observation  should  be  made.  If  Yn  Is  the  path 

.  A 

through  the  trellis  that  mislocates  the  edge  by  +n  points  and  Y-n  the  path  that 

mlslocates  the  edge  by  -n  points  then  In  general  C*  I1  C*  and 

YnY  Y-nY  A 

I  \Z'-  Y**  I  Ict^  i*  Hz"-  Y"*  I  |o^  .  Therefore  the  probability  that  Y  is  chosen 

■  i—  — n 1  y  — n  v  — n 

n  ,mn 

A 

over  Jf  will  not  be  the  same  as  the  probability  that  Y_n  Is  chosen  over  Y. 

This  lack  of  equality  between  Pr(n*a)  and  Pr(n“-a)  causes  the  random  variable 

n  to  have  a  nonzero  mean  value  which  introduces  a  bias  In  the  edge  location 

estimate.  Since  En  *  ZaPr(n-a)  and  since  Pr(n»a)  decreases  with  increasing 

a 

signal  to  noise  ratio  (SNR)  the  bias  will  decrease  with  increasing  SNR. 

Further  the  decision  boundaries  (and  hence  the  bias)  among  the  various 

2 

possible  paths  are  determined  by  the  variances  which  In  turn  are  deter¬ 
mined  by  the  possible  levels  a^,  and  the  degrading  function  h.  Thus  the 
probability  of  mlslocatlng  the  edge  by  n  points  (and  hence  the  bias)  Is 
Independent  of  the  sample  point  at  which  the  edge  occurred. 

If  the  SNR  is  high  enough  the  bias  can  probably  be  Ignored.  If  the  SNR 
Is  low  the  bias  can  be  calculated  by  the  procedure  mentioned  earlier  or  found 
experimentally  by  computer  simulation. 

Extension  to  Pulses 

In  measuring  the  size  of  an  object  two  edges  must  be  located.  If  the 
object  size  Is  sufficiently  large  or  the  signal  to  noise  ratio  Is  high  enough 


then  with  high  probability  the  minimum  cost  path  and  the  correct  path  through 
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the  trellis  will  coincide  somewhere  Inside  the  object. 

With  this  assumption  the  events  Ej  -  first  edge  mlslocated  by  n^  points 

and  “  second  edge  mlslocated  by  n^  points  are  Independent  [7].  Thus  If 

n^  Is  the  number  of  sample  points  difference  between  the  true  width  and  the 

estimated  width 

(12)  n  ■  n.  -  n„ 
w  2  i 

P(nw  -  a)  -  l  P(n2  -  a)  P(n^  -  a  -  8) 

0 

where  P(n2  ■  a)  and  P(n^  ■  a-  8)  can  be  calculated  from  (10).  The  sum  Is  over 
alt  possible  values  of  0  that  n2  can  assume  such  that  n^  -  n^  ■  a. 

If  the  noise  varies  with  the  signal  level  n^  may  not  have  zero  mean.  As 
In  the  previous  section  the  bias  can  probably  be  Ignored  if  the  SNR  is  high 
enough  or  calculated  theoretically  or  found  experimentally  If  necessary. 

Unknown  Levels 

In  the  previous  section  we  presented  formulas  for  the  probability  of  mls- 
locatlng  an  edge  by  |nj  points.  This  analysis  assumed  that  the  possible  levels 
aj  and  that  of  the  process  were  known.  In  this  section  the  probability  of 
error  Is  calculated  assuming  that  the  levels  are  not  known  a-priori. 

A  reasonable  course  of  action  In  the  case  of  unknown  levels  Is  to  obtain 
"training”  samples  of  the  gray  levels  characterizing  the  object  of  interest. 

A  A 

These  samples  can  be  processed  In  some  fashion  to  produce  estimates  a^  and  a^ 

A  A 

of  levels  aj  and  a2,  respectively.  Suppose  a^  ■»  a^  +  and  a2  "  a2  +  e2 

2 

where  and  e2  are  normal  random  variables  with  zero  mean  and  variance  Oj 

2  A  A 

and  0,,  respectively.  (If  aj  and  a2  are  large  sample  maximum  likelihood 
estimates  of  a^  and  a2  then  the  above  model  Is  accurate  [8],) 
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Let  the  correct  path  through  the  trellis  by  £  and  the  estimate  be  Y.  In 
this  case  the  correct  path  is  the  one  that  places  the  edge  in  the  correct 
position  even  though  the  levels  (or  height)  of  the  step  may  be  Incorrect.  The 
“signal  space"  diagram  of  Fig.  5  has  been  repeated  In  Fig.  6  with  the  addition 
of  a  new  point  0.  £  represents  the  output  sequence  corresponding  to  the  true 
edge  location  and  step  levels.  In  general.  £wlll  not  lie  in  the  plane 

A 

defined  by  Z,  Y,  Y.  However,  as  before,  only  the  projections  of  Y  -  and 
*  Y-Y 

Y  -  Z  where  Z  ■  0  +  N  on  -  will  have  any  effect  on  the  choice  between 

~  “  "  “  TTTW 

Y  and  Y. 

Now  Y-Y  *  Ah  *  (0,0, ...1,1,1, ...0,0)  A  ■  +(a.-a  ) 


n  l«s 


and 


Y-0  ■  h  *  (ej,e1»...,ej,e2,e2 . e^).  Thus  the  projection  T 


Y-Y 

of  Y-0  on  -  will  be  a  linear  combination  of  e.  and  e„.  Since  a  linear 

TF7TT  1  2 

combination  of  normal  r.v’s  Is  again  normal  the  only  Information  needed  to 

characterize  T**  is  its  variance  (it  will  have  zero  mean  since  and  were 

2 

assumed  to  have  zero  mean)  which  can  be  calculated  from  knowledge  of 

2  2 

,  o^,  n j ,  and  h» 

Again,  by  arguments  similar  to  those  in  the  previous  sections  the 

A 

probability  of  deciding  Y  over  Y  Is  the  probability  that  the  noise  along  the 

-  T 

Y-Y  axis  exceeds  -r-  “  T'or  equivalently  probability  that  the  sum  of  two  normal 

2  T 

2  2  1  n 

random  variables  with  variances  a  and  at  exceeds  — 


03)  -  QC—r—7 

2  Co2  +  a2) 


-)  and 


P  (E  )  -  2Q( - k  -n-  r  ) 

n  2K+°t> 
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This  result  can  be  extended  to  pulse  width  accuracy  as  before.  If 
•  Oj  -  n^  Is  the  number  of  sample  points  of  error  In  the  width  the  estimate 
then 

Pr[nw-o]  -  Z  Prfnj-B]  Pr(n2-a~3) 

6  T 

I  Q  A 

where  Pr(n.»3)  ■  Q( - r-  *— ).  Equation  (13)  shows  that  If  the  estimates  a, 

1  HoWt)  1 

■  *  2  C  2 

and  a2  are  such  that  c  < <  o  then  the  lack  of  a-prlorl  knowledge  about  the 
possible  levels  will  have  very  little  effect  on  the  probability  of  mlslocatlng 
an  edge. 

Simulation  Results 

A  pulse  of  width  thrlty  sample  points  was  generated  and  blurred  by  a 
linear  shift-invariant  system  with  a  Gaussian  shaped  Impulse  response  with  a 
standard  deviation  of  one  sample  point.  This  blurred  pulse  was  transformed  by 
y^  •  1,066  (log  1  -  1.5.)  +  .2  where  1^  Is  the  k—  sample  of  the  blurred 

pulse.  This  transformation  simulates  the  d-log  E  curve  of  film.  Noise  of 

1  r\ 

standard  deviation  .^Yk  J  was  added.  The  noisy  blurred  signal  was  then  pro¬ 
cessed  to  produce  an  estimate  of  the  pulse  width.  Several  different  SNR's 
were  used  with  one  hundred  width  estimates  obtained  at  each  SNR.  The  results 
are  shown  In  Figs.  6  and  7.  The  bias  In  the  width  estimate  can  be  seen  to 
decrease  with  Increasing  SNR  as  expected.  The  variance  of  the  estimate  can  be 
seen  to  decrease  with  Increasing  SNR  and  appears  to  asymptotically  approach 
the  Cramer-Rao  lower  bound. 
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DIGITAL  STRAIGHT  EDGES 
G.Y.  Tang  and  T.  S.  Huang 

Studies  of  the  digitization  of  a  strlght  line  have  been  made  by  Rosenfeld 
[1],  Freeman  [2],  Morse  [3l ,  Gaafar  [*♦] ,  Brons  [5].  A  set  of  rules  has  been 
established  to  govern  the  digitization  of  an  Ideal  straight  line  which  can  be 
described  by  a  first  order  mathematical  equation.  But,  In  the  real  world, 
most  edges  which  appear  to  be  straight  to  cur  eyes  do  not  fall  into  this  cate¬ 
gory.  On  the  contrary,  they  suffered  from  noise,  degradation,  blurring,  etc. 

In  this  report,  we  are  going  to  show  some  experimental  results  on  real  straight 
edges  and  to  compare  them  with  the  Ideal  case.  A  simple  testing  algorithm  is 
then  developed  to  determine  If  a  real  edge  Is  a  straight  one  or  not.  The 
application  of  the  testing  algorithm  can  be  found  in  various  areas  such  as 
syntactic  pattern  recognition,  character  recognition,  scene  analysis  and  ef¬ 
ficient  contour  coding,  etc.  Our  ultimate  goal  for  this  study  Is  to  use  sy- 
tactlc  method  to  aid  us  to  detect  edges  In  a  noisy  environment. 

In  the  following,  we  are  to  use  the  word  'real'  to  refer  to  whatever  Is 
on  pictures  obtained  by  practical  Imaging  devices.  The  word  'ideal'  refers  to 
the  Ideal  case. 

I.  PROPERTIES  OF  AN  IDEAL  STRAIGHT  LINE: 

An  Intensive  study  of  the  properties  that  a  digitization  of  a  straight 
line  should  have  has  been  reported  previously.  A  brief  summary  is  given  here. 

(A)  The  Chord  property,  proposed  by  Rosenfeld,  is  a  necessary  and  suffi¬ 
cient  condition  for  a  digital  arc  to  be  straight. 

(B)  The  digitization  of  any  real  curve  can  be  expressed  In  terms  of  a 
chain  code  [2],  The  chain  code  of  a  straight  line  should  obey  the  following 
rules. 

1)  There  are  at  most  two  different  code  elements  differing  by  1,  modulo  8. 
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2)  One  of  these  two  codes  occurs  singly. 

3)  Successive  occurrences  of  the  single  element  are  as  uniformly  spaced 
as  possible. 

C)  For  an  Ideal  straight  line  with  rational  slope,  the  corresponding  chain 
code  Is  a  repetition  of  a  certain  period. 

D)  We  can  always  find  a  rational  slope  to  approximate  a  real  slope  so  that 
their  chain  codes  are  as  close  to  each  other  as  we  want. 

E)  An  ad  hoc  testing  algorithm  has  been  proposed  [3]  to  see  If  a  given 
chain  code  may  come  from  a  straight  line. 

II.  EXPERIMENT 

The  nice  properties  listed  In  the  foregoing  section  are  based  on  the 
assumption  that  the  lines  are  Ideally  straight,  the  digitization  Is  noise-free 
and  the  digitization  scheme  Is  well  defined  [l],  [2],  In  the  real  world,  where 
a  model  with  such  a  high  order  of  Idealism  can  hardly  find  many  applications, 
minor  randomness  corrupts  the  structures  severely.  The  Intention  of  our  ex¬ 
periment  Is  to  observe  the  code  structures  of  real  edges  and  to  compare  them 
to  the  Ideal  case.  Aiso  we  outline  a  new  testing  algorithm  from  our  experi¬ 
mental  results. 

A  set  of  12  pictures  has  been  taken  by  a  standard  Nikon  F-2  camera  with 
35  mm  film.  A  flying-spot  scanner  Is  used  to  discretize  the  picture  Into  a 
square  matrix  of  size  256  by  256.  The  grey  levels  range  from  0  to  255*  A 
simple  contour  follower  and  a  chain  code  encoder  Is  used  to  obtain  the  chain 
codes  on  each  picture.  Only  straight  edges  are  on  these  12  pictures.  There 
are  totally  30  straight  edges. 

A  contour  follower  Is  designed  so  that  the  output  of  the  contour  follower 
Is  the  chain  codes  of  the  edge  which  It  follows.  The  input  to  the  contour 
follower  Is  a  window  of  size  3x3.  We  assume  that  the  center  of  the  window  Is 
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on  the  contour  and  that  we  know  which  of  Its  eight  neighbors  Is  also  on  the 
contour.  The  task  of  the  contour  follower  Is  to  find  a  point  which  Is  supposed 
to  be  the  next  contour  point  among  the  eight  neighbors  around  the  center.  Then 
the  window  moves  to  and  centers  at  that  point.  The  direction  of  move  Is  en¬ 
coded  by  Freeman1 s  method.  Repeating  the  same  procedure,  we  can  obtain  the 
whole  contour.  The  way  to  locate  the  next  point  from  the  eight  neighbors  of 
the  center  Is  simply  looking  at  the  eight  neighbors  In  counter-clockwise  sense 
starting  with  the  neighbor  already  known  on  the  contour.  A  thresholding  tech¬ 
nique  Is  used  to  discriminate  the  two  regions  defining  the  contour. 

Table  1  shows  which  ideal  properties  are  violated  by  chains  obtained  from 
real  pictures. 

The  Ideal  edge  Is  not  what  one  finds  In  the  Images  produced  by  real-life 
imaging  devices.  There  are  several  factors  that  degrade  the  edges  that  are 
actually  found.  Two  predominant  ones  are  blurring,  or  defocusing,  and  Ir¬ 
regularities  of  the  surface  structure  of  the  object.  Besides,  for  the  case  of 
straight  edges,  a  slight  concavity  or  convexity  can  be  considered  as  another 
kind  of  disturbance.  The  net  effect  of  these  disturbances  of  Freeman  chain 
code  can  be  summarized  as: 

1)  A  third  code  element  Is  introduced  If  there  Is  a  missing  or  spurious 
lattice  point.  This  third  code  element  will  occur  together  with  the  code 
element  which  occurs  singly.  Figure  I  illustrates  this. 

2)  A  third  run  length  Is  Introduced  If  the  missing  or  spurious  lattice 
point  occurred  at  the  beginning  or  at  the  end  of  each  run.  Fig.  2  Illustrates 
this. 

3)  Other  types  of  errors  occur  which  cannot  be  characterized  easily. 

Fig.  3  shows  us  the  test  pictures. 
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III.  STRAIGHT  LINE  RECOGNITION 

To  recognize  a  real  chain  code  as  a  straight  line  Is  not  straightforward 
since  we  have  no  rigorous  definition  of  what  is  a  real  straight  edge.  Here  we 
propose  to  employ  a  heuristic  method  to  solve  the  problem  of  the  recognition 
of  digital  straight  lines. 

The  basic  strategy  of  our  approach  Is  rather  simple.  We  ar  attempting  to 
enclose  a  digital  straight  line  by  a  straight  band.  If  such  a  band  with 
reasonable  bandwidth  could  be  found  for  a  given  chain  code,  then  we  claim  that 
that  chain  code  corresponds  to  a  straight  line.  Mathematically,  a  straight 
band  is  defined  by  an  equality 

|  Y  -  mX|j<_k 

where  m  Is  the  slope  of  the  straight  band  and  k  is  related  to  the  bandwidth  of 
the  straight  band.  Two  parameters,  m  and  k„  are  therefore  to  be  determined 
from  a  given  chain  code.  We  proceed  as  follows: 

Step  1  ■  Check  If  there  are  only  two  code  elements;  one  occurs  singly. 

Yes;  go  to  Step  3 
No;  go  to  Step  2. 

Step  2  ■  Correct  all  possible  first  kind  of  errors  mentioned  In  the  pre¬ 
vious  section.  Then  check  If  there  are  only  two  code  elements  left;  one 
occurs  singly,  if  It  Is  yes,  then  go  to  Step  3;  otherwise  do  the  following: 

Calculate  the  percentage  of  the  3rd,  4th  ...  etc.  code  elements  to  the 
total  length  of  the  chain  code  and  the  percentage  of  the  second  code  elements 
which  should  occur  singly  in  Ideal  case  but  occurs  in  runs  in  this  case  to  the 
total  number  of  the  second  code  elements.  Then  compare  these  two  percentages 
to  two  preset  thresholds.  Once  they  are  smaller  than  the  thresholds,  go  to 
Step  3;  otherwise  the  chain  code  Is  rejected  as  a  straight  line. 
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Step  3  *  m  Is  the  ratio  of  the  element  which  occurs  singly  to  the  total 
number  of  1st  and  2nd  elements. 

k  Is  a  subjective  quantity.  It  relates  to  the  noise  level  of  the  picture, 
the  sharpness  of  the  edge  and  the  resolution  of  the  picture. 

It  seems  that  so  far  we  have  taken  few  advantages  of  the  nice  structural 
regularities  existed  among  ideal  straight  lines  to  aid  us  in  developing  test* 
ing  algorithm  for  the  non-ideal  case.  A  rather  straightforward  heuristic 
method  has  been  devised  to  further  reduce  the  number  of  testing  points.  It 
Is  noted  that  the  straight  band  is  a  geometrically' convex  set.  So,  only  the 
two  ends  of  each  run  on  the  chain  codes  need  to  be  tested.  Furthermore,  if  we 
apply  a  simple  operation  which  is  to  replace  each  run  by  its  run  length  and  to 
delete  the  code  elements  occurred  singly,  then  the  output  of  the  operation 
should  obey  the  rule:  only  two  code  elements;  one  occurs  singly,  for  the 
Ideal  case.  For  the  real  chain  code,  we  may  apply  the  same  operation  itera¬ 
tively  to  It  until  the  output  fails  to  follow  the  rule.  We  then  use  the 
final  one  as  a  clue  to  break  the  initial  input  chain  Into  several  pieces  oF 
line  segments.  Only  the  ends  of  each  line  segment  are  to  be  tested.  An 
example  illustrating  this  follows: 

A  chain  code  from  picture  No.  7  in  our  experiment  is 

A 

010101010101010101010100101010101010101 

A  A 

010101010100101010101010101010010101010 

A  A 

101010101019100101010101010101010010101 

A 

01010101010101001010101010  1. 
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The  first  output  of  our  operation  Is: 

A  A  A 

111111II11I2111M11111112M11I1112 

A  A  A 

I  I  1  1  l  l  1  1  1  1  2  1  1  '  1  1  1  1  1  2  1  1  1  1  1  1  1  1  1  2  1  1  1  1  1. 

The  next  output  Is: 

II  12  8  10  8  9  5  .  It  falls  to  follow  the  rule.  So 

we  trace  back  and  locate  the  potential  testing  points  (wlthA  )  on  each  gen¬ 
eration.  The  coordinates  of  the  testing  points  with  respect  to  the  Inltlum 
are  then  (22,11),  (49,24) ,  (68,33),  (90,44) ,  (110,53),  (131  ,63). 

Table  2  shows  the  result  of  applying  the  foregoing  method  to  test  17 
chain  codes  obtained  In  our  experiment.  The  "MIN  THRESHOLD"  means  the  small¬ 
est  k  value  we  have  to  choose  In  order  to  report  the  edges  straight.  Also  we 
applied  Morses  algorithm  to  test  these  real  chain  code.  It  does  not  work  as  well. 
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)•  Only  two  code  elements 

2.  One  of  the  two  code  elements  occurs  singly 

3.  Only  two  run  lengths 

A.  Two  run  lengths  are  consecutive  integers 
5*  ®ne  the  two  run  lengths  occurs  singly. 

6.  The  runs  of  the  run  lengths  can  be  of  at  most  two 
run  lengths:  one  of  which  occurs  singly 


Table  1.  Comparison  of  real  cases  and  ideal  case 
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Ltne  No.  length 


3-1 

147 

3-2 

49 

3-3 

54 

3-* 

58 

4-1 

82 

4-2 

*7 

4-3 

42 

4-4 

127 

5 

238 

6 

237 

7 

130 

8 

209 

9-1 

85 

9-2 

245 

10-  r 

147 

10-2 

40 

11-2 

181 

Min.  Threshold  MORSE 


.795 

NO 

.755 

NO 

.833 

NO 

.724 

NO 

1.329 

NO 

.404 

YES 

1.524 

NO 

1.055 

NO 

1.345 

NO 

1.456 

NO 

.435 

YES 

1.004 

NO 

1.17 

NO 

.942 

NO 

1.258 

NO 

.900 

NO 

1.127 

NO 

Table  2  The  result  of  the  proposed  method  In  comparison 
with  Morse  algorithm. 
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perfect  case  : 


X-X-X-X-X-X-X-X-X 


X-X-X-X-X-X-X-X-X 


00000000100000000 


e  missing  point: 


X-X-X-X-X-X-X-X 


X-X-X  X-X-X-X 


0071000100000000 


two  missing  point: 


X-X-X-X-X-X-X-X 


X-X-X 


X-X-X 


X  -  X 


0  0  701001000000 


a  spurious  point: 


X-X-X-X-X-X-X-X 


X-X-X  X  X-X-X-X 


001  7000100000000 


two  spurious  point: 


X  -  X 


X-X-X-X-X-X-X-X 


X-X-X  X  X  X-X-X 


0010700100000000 


Fig.  1.  Hissing  or  spurious  points  occurred  In  runs 


perfect  one 


3* 


X  -  X  -  X  -  X 

X  -  X  -  X  -  X 

X  -  X  -  X  -  X 
00010001000 

X  -  X  -  X  -  X 

a  spurious  point: 

X-X-X-X-X 

X  -  X  -  X  -  X 
00100001000 
a  missing  point: 

X  -  X  -  X  -  X 

X  -  X  -  X 

X-X-X-X-X 

00001001000 

Fig.  2  Missing  or  spurious  points  occured  at  the 
ends  of  a  run 


.  5 
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IMAGE  STRUCTURE:  SYNTACTIC  SCENE  ANALYSIS 
K.  S.  Fu  and  R.  Y.  LI 

The  essential  problem  In  our  research  on  the  syntactic  pattern  recognition 
of  highways  and  rivers  is  really  to  find  a  grammar  that  will  describe  well 
these  classes  of  Interest.  If  the  physical  shape  of  the  class  under  consider¬ 
ation  is  completely  known  and  fixed,  like  printed  English  character;  we  can 
Immediately  write  down  the  syntactic  rules  to  describe  Its  structure.  Since 
this  Is  not  quite  true  In  our  case,  the  construction  of  grammatical  rules  has 
to  be  based  on  Informations  obtained  from  a  set  of  sample  patterns  known  to 
come  from  that  class.  Hopefully,  this  set  of  Inferred  rules  should  be  able  to 
describe  and  predict  other  sample  patterns  which  are  of  the  similar  nature  as 
the  original  training  samples  and  presumably  In  the  same  class.  A  basic  ap¬ 
proach  of  grammatical  Inferrence  problem  is  to  construct  a  grammar  by  Identify¬ 
ing  the  syntactic  structures  of  the  known  string  and  any  possible  recursiveness 
that  might  happen.  There  are  three  steps: 

(1)  Try  to  discover  the  syntactic  structure  of  the  given  string  by  looking  for 
repetition  and  dependent  relationships. 

(2)  Decide  what  sublanguages  make  up  the  language  and  generate  non-terminals 
for  each  sublanguages. 

(3)  Combine  equivalent  nonterminals  which  have  almost  the  same  sublanguage  and 
determine  the  appropriate  relationships  among  sublanguages. 

One  practical  method  to  learn  the  syntactic  structures  of  the  given  pictures 
is  to  use  a  semantic  teacher  to  learn  the  meaningful  nonterminals  one  level  at 
a  time  [1],  To  start  the  Inferrence  process,  we  first  find  the  types  of  term¬ 
inals  or  primitives  that  will  fit  the  subparts  of  the  picture  pattern  for  a 
given  window  size.  After  this  initial  extraction  process,  we  have  to  decide 
the  most  probable  combinations  of  primitives  which  occur  as  neighbors  of  each 
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other  in  the  set  of  observed  training  samples.  These  combinations  are  then 
applied  to  the  training  data  set  to  test  their  recognition  effectiveness.  When 
the  results  appear  tc  be  satisfactory  after  some  additions  and  deletions  of  the 
combination  rules,  we  can  choose  this  set  of  rules  to  represent  the  training 
samples.  The  appropriate  grammar  can  then  be  formulated  by  using  these  rules. 
In  our  present  case,  we  chocse  the  tree  grammar  because  of  its  easiness  to 
describe  these  rules.  A  tree  recognition  program  based  on  this  tree  grammar 
can  then  be  used  to  recognize  the  training  data  set  of  Lafayette  and  a  test 
data  set,  that  of  Grand  Rapids,  Michigan,  Figures  (1),  (2),  and  (3)  contain 
the  results  from  Lafayette  experiment.  Other  preliminary  results  from  Grand 
ftaplds  can  be  found  In  reference  [3l. 
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Figure  (1)  Statistically-classified  results  of  Lafayette 


*♦0 


•  •  »r 

*■*  *•  *•  « 


M  •  •  « 

•HIINI*  «* 


V 

■  • 


.•*  :  : 


-■  •  .  j 


«•  .  *. 
M«  •  ft 

•  I  Ml 


N  M  w 

.  -  . 


•  H  ft 

J  *• 

!2  ••  •• 


p’  p’ 


Figure  (2)  Skeltons  from  Lafayette  samples 
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SYNTACTIC  SCENE  SEGMENTATION 
K.  S.  Fu  3rd  J.  Keng 

A  syntactic  approach  to  scene  segmentation  has  been  investigated  which 
Involves  two  levels  of  processing.  The  first  level,  referred  to  as  the 
transformation  process,  consists  of  five  steps  referred  to  as  (1)  threshold 
finding,  (2)  horizontal  processing,  (3)  vertical  processing,  (4)  logic 
Integrating,  and  (5)  line  smoothing.  The  second  level,  which  is  the  actual 
syntactic  analysis,  requires  Inference  of  a  tree  grammar  to  describe  the 
boundaries  of  homogeneous  regions.  The  tree  grammar  is  then  Implemented  in  a 
parser  which  traces  the  region  boundaries. 

The  approach  has  been  implemented  and  initial  experiments  on  multi- 
spectral  remote  sensing  imagery  are  being  conducted.  Further  detail  and 
experimental  results  will  follow. 

Working  on  the  problem  of  picture  segmentation  through  a  syntactic 
approach,  we  feel  that  the  evaluation  of  the  earth  resources  Is  very  useful. 

So,  the  multispectral  remotely  sensed  picture  Is  chosen  as  data. 

There  are  two  levels  of  processing  for  the  syntactic  picture  segmenta¬ 
tion,  first,  the  transformation  process  and,  second,  the  tree  grammar  analysis. 
The  first  level  consists  of  five  sub-processes,  threshold  finding,  horizontal 
processing,  vertical  processing,  logic  Integrating,  line  smoothing.  The 
process  of  tree  grammar  analysis  utilizes  the  corresponding  parser  from  the 
Inferred  tree  grammars  to  process  the  transformed  picture.  Then  a  picture  Is 
segmented. 
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Syntactic  Picture  Segmentation 
a.  First  level  -  transformation  processing 

(1)  Threshold  finding 

From  a  digitized  picture,  training  regions  are  located  from  every 
homogeneous  parts  of  the  picture,  then  the  differences  between  two  of 
these  means  of  grey  levels  of  the  training  regions  are  calculated.  If 
they  are  close  to  zero  we  neglect  them,  then  the  minimum  of  them  Is 
selected  as  a  threshold. 

(2)  Horizontal  processing 

The  model  of  the  multisectral  images  Is  defined  on  the  Euclidean 

n-d (mens Iona  1  space  En  which  Is  a  space  having  n  coordinates 

X|,X. . xn>  The  number  n  here  represents  the  number  of  channels  to 

be  chosen.  A  point  of  the  space  Is  by  definition  an  ordered  n-tuple 

(x, ,x_ , • • • ,x  ).  The  distance  between  two  points  is  defined  as 
i  c  n 

Euclidean  distance.  The  operation  procedure  Is  to  compare  the  grey 
levels  of  (1,1)  and  (1,2).  If  the  distance  is  smaller  than  threshold, 
then  set  zeroes  to  (1,1)  and  (1,2).  Then,  (1,1)  and  (1,3)  are  compared, 
if  the  distance  Is  greater  than  threshold.  A  one  is  put  to  (1,3)  and 
the  same  operations  start  from  (1,A).  The  operations  on  second  row 
follow  the  same  pattern. 

(3)  Vertical  processing 

The  operations  are  the  same  as  horizontal  processing  except 
It  goes  vertically  instead  of  horizontally. 

(k)  Logic  Integrating 

The  logic  variable  of  horizontal  processing  Is  named  as  H  and  V 
for  the  vertical  processing.  .The  ioglc  integrating  process  achieves 
the  Integration  through  Boolean  algebra  V+H. 
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(5)  Line  smoothing 

The  line  smoothing  algorithm  connects  the  discontinuity  of  the 
lines.  As  for  this  part,  lots  of  line  smoothing  algorithms  can  be 
devised. 

b.  Second  level  -  tree  grammar  analysis 

A  tree  grammar  Is  Inferred  to  describe  the  boundaries  of  the 
homogeneous  regions.  The  tree  grammar  Is  used  to  trace  the  boundaries 
and  reject  the  unnecessary  boundary  parts.  Finally  a  segmented  picture 
Is  received  from  the  (two  level)  syntactic  method. 

The  scheme  of  syntactic  picture  segmentation  has  been  Implemented  and 
the  experiments  have  been  conducted  on  the  multispectral  remotely  sensed 
pictures  of  an  area  In  the  State  of  Indiana.  The  pictures  were  taken  on 
August  13,  1971  and  stored  In  the  computer  IBM  360/67  of  the  Laboratory  for 
Applications  of  Remote  Sensing  In  West  Lafayette,  Indiana. 

The  result  of  a  picture  96x96  area  Is  shown  In  Fig.  1.  For  the  purpose 
of  comparing  processing  time  and  accuracy,  a  result  of  statistical  segmenta¬ 
tion  by  clustering  Is  also  provided  In  Fig.  2.  The  computer  processing  time 
on  IBM  360/67  for  the  same  area  of  the  picture,  the  syntactic  picture  seg¬ 
mentation  takes  about  36  seconds  and  the  clustering  technique  requires  180 
seconds. 

The  area  of  column  (1 05”1 32 )  and  row  (444-450)  In  Fig..  I,  shows  lots  of 
boundary  points.  From  a  survey  of  the  ground  truth  (Fig.  3),  It  points  out 
the  reason.  Because  this  area  Is  pasture  and  bare  soli  compound  area.  So 
lots  of  boundary  between  these  two  objects  of  the  compound  area  are  located. 
A  classified  result  of  the  same  area  Is  also  provided  In  Fig.  4.  In  the 
region  column  (90-100),  row  (480-488)  In  the  syntactic  segmentation  result 
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shows  a  segment  which  corresponds  to  the  same  location  in  the  classified 
result.  But  the  clustering  result  can  not  segment  it  well. 
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Figure  1  The  syntactic  picture  segmentation  result 
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Figure  2  The  statistical  picture  segmentation  Clustering) 
on  same  area  as  Fig.  1 
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Figure  4  The  classified  result  of  the  area  of  Fig.  1 


STATISTICAL  DEPENDENCY  MODELS  OF  CONTEXT 
P.  H.  Swain  and  T.  S.  Yu 

In  the  previous  quarterly  report,  we  described  a  model  for  Incorporating 
context  In  the  Image  analysis  process  through  compound  decision  theory. 
Briefly,  under  a  fairly  stringent  set  of  assumptions  a  Bayesian  strategy  Is 
employed  which  classifies  a  point  Into  one  of  a  candidate  set  of  classes  based 
on  the  mult  I  spectral  data  from  the  point  itself  and  the  data  from  either  the 
neighboring  four  or  neighboring  eight  points.  Figure  1  shows  the  results  of 
applying  this  approach  for  classifying  a  small  set  of  LANDSAT  mul tispectral 
scanner  data,  A  block  of  Imagery  128x128  pixels  was  classified  using  the 
‘Simple11  rule  (no  context),  with  the  ^-neighbor  rule,  and  with  the  8-neighbor 
rule.  Samples  of  900  pixels  (30x30)  were  selected  exhibiting  a  range  of 
polntwlse  accurfes  and  the  corresponding  accuracies  obtained  using  the  context 
Incorporating  rules  were  tabulated  together  with  the  corresponding  classifica¬ 
tion  times.  The  results  of  the  experiment  show,  as  expected,  the  classifier 
using  context  Is  consistently  better  than  the  classifier  which  makes  each 
decision  based  on  data  from  the  individual  points.  Furthermore,  the  "8- 
nelghbor  classifier"  Is  consistently  better  than  the  "4-neighbor  classifier". 
However,  the  price  paid  in  terms  of  computation  time  is  substantial.  To 
Justify  general  use  of  this  approach  we  shall  have  to  demonstrate  (a)  Its 
performance  potential  over  a  sufficiently  broad  range  of  analysis  problems, 
and  (b)  a  means  of  Implementation  (special  purpose  hardware,  software,  or  a 
combination  thereof)  which  is  efficient  enough  to  provide  results  on  a  cost- 
effective  basis. 
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IMAGE  NOISE  REDUCTION 
M.  Y.  Yoo  and  I.  S.  Huang 

I.  INTRODUCTION 

There  have  been  many  attempts  to  enhance  images  degraded  by  detector  noise 
and  Imperfection  of  the  imaging  system  [1-4].  In  most  cases  we  encounter  two 
contradicting  requirements:  reducing  the  noise  as  much  as  possible,  and  retain¬ 
ing  edge  sharpness.  The  only  reasonable  answer  !s  a  compromise  between  the 
two.  We  will  approach  this  problem  using  the  synthetic  highs  technique  where 
we  can  treat  the  smoothly  varying  part  and  the  sharply  changing  boundaries 
separately  and  we  may  enjoy  some  freedom  in  compromising  the  two  situations 
t5].  The  performances  of  the  system  will  be  compared  with  several  available 
heuristic  approaches  [6],  [7]. 

I I .  IMAGE  ENHANCEMENT  SYSTEMS 

2. 1  Synthetic  Highs  Technique 

This  technique  was  proposed  by  Schreiber  [8]  and  was  used  in  two  dimen¬ 
sional  contour  coding  for  data  compression  by  Graham  [9] .  The  basic  idea  of 
this  technique  is  to  decompose  the  image  into  two  major  elements  (slowly 
varying  "lows"  and  synthetic  highs"  which  sre  mostly  boundaries  and  textured 
parts)  such  that  the  recombination  of  the  two  elements  results  in  the  original 
again. 

The  technique  was  described  in  detail  in  previous  reports  and  In  Graham 
[9],  and  we  are  not  going  to  repeat  that  here  again.  To  use  the  technique  for 
noise  reduction,  we  need  a  noise  reduction  filter  between  the  edge  detector 
and  the  reconstruction  filter.  So  we  need  an  acceptable  boundary  detector  for 
noisy  Images.  Tang's  [10]  or  the  following  simple  edge  detector  may  be 
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Noise  Reduction 
Filter 


Reconstruction 
F! iter 


Reconstructed 

Picture 


Figure  1  Synthetic  highs  technique 

First  we  calculate  the  grey  level  differences  in  four  different  dlrec 
tions;  horizontal,  vertical,  diagonal  1  (**5°)»  and  diagonal  2  (135°)  as 
fol lows: 


G(x-Ax.y)  -  G (x+Ax,y) I 
e  G(x-Ax,y) +  G(x+Ax|yj 


dlag  2 


|G(x-Ax,y-Ay)  -  G(x+Ax,y+Ay| 
G(x-Ax,y-Ay)  +  G (x+Ax,y+Ay) 


Whenever  one  of  the  following  conditions  holds  we  decide  Gfx,y)  is  on  the 
boundary. 

1)  Horz  0j  and  Vert  <  ©2 

11)  Horz  <  0?  and  Vert  _>  9^ 

III)  diag  1  0j  and  diag  2  <  0^ 

Iv)  diag  2  <  02  and  diag  2  >.  0j 

where  ®2  ®2  —  are  Preass '9nec*  threshold  values.  Typical 
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2.2  Median  Filter 

We  use  an  n  x n  (n  is  odd)  window  to  scan  the  noisy  image  and  replace  the 
grey  level  of  the  center  by  the  median  grey  level  within  the  window.  Pratt 
[6]  used  one-dimensional  window  to  remove  impulse  noise.  Resolution  of  the 
filtered  Image  is  highly  dependent  upon  the  size  of  the  window  and  appropriate 
size  should  be  chosen  to  retain  reasonable  boundaries.  3x3  and  5x5  windows 
were  used  for  our  experiment. 

2.3  Variable  Width  Filter 

First  we  take  the  gradient  of  the  noisy  image  and  divide  the  absolute 
gradient  level  into  four  different  Intervals.  We  assign  the  smooching  filter 
of  an  appropriate  size  to  each  interval  such  that  the  lower  the  absolute 
gradient  level  Is,  the  larger  the  duration  of  the fi  Iters  impulse  response  Is. 
One  dimensional  Gaussian  and  median fi  1  ters  can  be  used  in  the  horizontal  and 
the  vertical  directions  sequentially. 

The  basic  idea  of  this  approach  Is  that  we  retain  more  boundaries  where 
the  absolute  gradient  level  is  high  by  using  narrow  smoothing  filters  and 
heavily  smooth  out  slowly  varying  part  by  wide  smoothing  filters.  Uniform  or 
non-uniform  subdivision  of  gradient  levels  may  be  used  depending  upon  the 
distribution  of  gradient  levels.  The  sizes  of  the  filters  used  areO,  3,  5,  7. 

2.4  Noise  Cheating  Technique 

This  technique  Is  a  combination  of  two  averagings  with  different  window 
sizes.  We  average  the  noisy  image  using  an  mxm  window  and  average  the  noisy 
picture  again  by  an  n x n  (n  >  m)  window.  In  the  original  paper  [7]  the 
authors  quantize  grey  levels  of  "severely"  averaged  images  using  a  quantum 
step  that  is  at  least  four  times  the  standard  deviation  of  the  averaged 
picture.  But  if  the  standard  deviation  of  the  averaged  picture  is  large 
enough,  the  quantizing  process  wipes  out  everything  and  this  really  happened 


In  our  case.  So  we  may  have  to  skip  the  quantizing  step  depending  upon  the 
size  of  the  standard  deviation.  After  averaging  the  noisy  Image  with  two 
different  sizes  of  windows  we  replace  each  grey  level  In  the  "lightly"  averaged 
Image  by  the  closest  grey  level  among  the  corresponding  eight  surrounding 
Image  points  In  the  "severely"  averaged  Image. 

The  Idea  Is  that  we  retain  the  resolution  of  the  "lightly"  averaged 
picture,  while  we  enjoy  the  reduction  of  noise  level  of  the  "severely"  averaged 
Image.  The  combining  scheme  is  a  kind  of  discrete  maximum  likelihood  approach. 
m«2,  n«3  were  used  for  the  experiment. 

III.  EXPERIMENTAL  RESULTS  AND  CONCLUSION 

White  Gaussian  noise  was  added  to  generate  a  10  dB  (variance  of  signal/ 
variance  of  noise)  noisy  picture  of  size  256x256.  The  bandwidth  of  the  low 
pass  Gaussian  filter  in  synthetic  highs  system  is  0.116.  The  original  noise¬ 
less  picture  is  shown  In  Fig.  2  and  the  10  dB  noisy  picture  and  "lows"  of  the 
noisy  In  Fig.  3  and  Fig.  4,  respectively.  Figure  5  shows  the  noise  reduction 
filter  used  In  synthetic  highs  system.  But  this  filter  should  be  extended  by 
a  3x3  window  so  as  to  pass  both  the  positive  and  the  negative  parts  of  boun¬ 
daries  detected  by  Laplacian  or  gradient  edge  detector,  otherwise  we  lose 
resolution  significantly.  Reconstructed  pictures  are  given  in  Fig.  6  and 
Fig.  7.  Outputs  of  the  two-dimensional  median  filter  are  shown  In  Fig.  8  and 
Fig.  9.  In  5x5  median  filtered  picture  we  lost  resolution  quite  a  bit. 

Median  filters  seem  to  be  very  effective  for  removing  impulse  errors  (all  dark 
spots  have  gone  away  in  both  outputs),  but  still  retain  Gaussian  noises  at  an 
unpleasant  level.  Variable  width  Gaussian  filters  are  truncated  at  2  x 
standard  deviation. 


A  one-dimensional  filter  was  used  sequentially  In  horizontal  and  vertical 
directions  and  when  we  apply  the  filter  in  the  vertical  direction  we  use  data 


already  averaged  In  the  horizontal  direction.  This  may  cause  more  reduction 
of  noise  but  resolution  Is  lost  also.  Actually  we  didn't  smooth  out  at  the 
top  Interval  where  the  absolute  gradient  Is  highest-.  (Filter  size  was  specified 
by  0  In  section  2.3.)  Comparing  with  the  variable  width  Gaussian  filter, 

the  performance  of  the  variable  width  median  filter  Is  very  poor.  In  the  noise 
cheating  technique  we  first  averaged  the  noisy  picture  by  2x2  window  and  re¬ 
placed  the  grey  levels  at  4  picture  points  in  the  window  by  the  averaged  level. 

We  did  the  same  thing  with  3x3  window  and  el  Imlnated  Isolated  picture  blocks 
(Note:  3x3  window  will  have  the  same  grey  level  after  averaging)  by  simply 
replacing  the  center  grey  level  by  the  surrounding  grey  levels.  Since  the 
whole  block  (2x2, 3x3)  has  the  same  grey  level,  the  output  of  the  noise  cheating 
technique  has  lots  of  square  blocks.  The  performance  of  this  technique  seems 
the  worst.  The  two-dimensional  median  filter  is  most  economical  and  easiest 
to  apply  but  the  synthetic  highs  system  gives  the  best  result  although  it 
Is  most  costly. 

Some  of  the  techniques  can  be  modified  for  better  performance.  For 
example,  we  may  use  the  noise  reduction  filter  used  in  the  synthetic  highs 
system  for  variable  width  Gaussian  or  median  filter  and  we  may  use  moving  3x3 
overlapping  window  replace  the  center  grey  level  by  the  averaged  level  rather 
than  assigning  the  same  grey  level  to  the  whole  block.  Finally,  we  emphasize 
that  the  performance  of  noise  reduction  techniques  depends  upon  the  type  of 
noises  involved. 

REFERENCES 

[1]  T.  S.  Huang  and  et.  al.,  "Image  Processing,"  Proc.  IEEE,  Vol.  59,  No.  11, 
November  1971.  pp.  1586-1609. 

— — ,  IEEE  Proc.,  Vol.  60,  Special  Issue  on  Digital  Image  Processing, 

July  1972. 


[2] 


57 


[3l  W.  K,  Pratt,  "Generalized  Wiener  Filtering  Computation  Techniques,"  IEEE 
Trans,  on  Comp.,  Vol.  C-21,  No.  7,  July  1972,  pp,  636-641. 

[4]  A.  Habib! ,  "Two-Dimensional  Bayesian  Estimation  of  Images,"  Proc.  IEEE, 
Vol.  60,  No.  7,  July  1972,  pp.  878-883. 

[5]  M.  Y.  Yoo  and  T.  S.  Huang,  "Noise  Reduction  In  Photographic  Images," 
RADC-TR-75-307,  PP.  50-54. 

[6]  W.  K.  Pratt,  "Median  Filtering,"  USCIPI  Report,  pp.  116-122. 

[7]  H.  J.  Zwelg,  E.  B.  Barrett,  and  P.  C.  Hu,  "Noise  Cheating  Image  Enhance¬ 
ment,"  JOSA,  Vol.  65,  No.  11,  November  1975,  pp.  1347-1353. 

[8]  W.  F.  Schrelber,  "The  Mathematical  Foundation  of  the  Synthetic  Highs 
Systems,"  Research  Lab  of  Electronics,  M.I.T.,  OPR,  Vol.  68,  January 
1963,  pp.  140. 

[9]  0.  N.  Graham,  "Image  Transmission  by  Two-Dimensional  Contour  Coding," 
Proc.  IEEE,  Vol.  55,  No.  3,  March  I967. 

[10]  G.  Y.  Tang  and  T.  S.  Huang,  "Edge  Extraction  Technique,"  RADC-TR-75-202, 
August  1975. 


.tsj*«iv  .'-V,  .-.  „■>  jvgvVVV 


.vv-. 


60 


Figure  6  Reconstructed  picture 
Leplecien  synthetic 
highs  used 


Figure  7  Reconstructed  picture 
Gradient  synthetic 
highs  used 
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Figure  8  Reconstructed  picture 
3x3  median  filter  used 


Figure  9  Reconstructed  picture 
5x5  median  filter  used 


Figure  10  Reconstructed  picture 
Variable  width  Gaussian 
uni fora  quantization  used 


Figure  11  Reconstructed  picture 
Variable  width  Gaussian 
non-uniform  quantization 
used 
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Figure  12  Reconstructed  picture 
Variable  width  Median 
filter  uniform 
quantization  used 


Figure  13  Reconstructed  picture 
Variable  width  median 
filter  non-uniform 
quantization  used 
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Ftgur*  U  flecon«tructed  picture 

Noise  cheating  technique 
2x2  and  3x3  windows  used 
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TWO-DIMENSIONAL  DIGITAL  RECURSIVE  FILTERS  AND  THEIR 
APPLICATION  TO  IMAGE  PROCESSING 

Brian  O'Connor  and  T.S.  Huang 

IWo  dimensional  digital  recursive  filters  have  the  potential  of  saving 
computer  time  and  storage  In  the  processing  of  large  two  dimensional  arrays  such 
as  Images.  In  our  research  we  are  concerned  with  their  design  and  their  appli¬ 
cation  to  images.  The  desired  processing  to  be  performed  on  an  Image  can  gen¬ 
erally  be  described  In  the  frequency  domain.  In  two  dimensional  filter  design 
the  filter  coefficients  are  found  which  best  approximate  this  desired  frequency 
response,  whMe  at  the  same  time  guaranteeing  stability. 

Several  possible  design  approaches  exist.  One  relies  on  nonlinear  optimi¬ 
zation  techniques  to  minimize  the  l  norm  of  the  difference  between  the  desired 
magnitude  or  group  delay  response  and  the  filter  response  [1,2],  These  tech¬ 
niques  can  be  modified  to  guarantee  stability.  However,  numerical  problems 
arise  with  nonlinear  optimization  techniques,  namely  very  large  amounts  of 
computation,  sensitivity  to  starting  points,  and  the  possibility  of  converging 
to  a  local  rather  than  global  minimum.  In  addition.  It  Is  necessary  to  check 
stability  of  the  designed  filter  at  each  iteration  of  the  algorithm.  However, 
a  recently  proposed  method  has  eliminated  need  for  the  stability  tests  [3]. 
Another  approach  Is  based  on  an  algorithm  which  allows  the  designer  to  specify 
an  arbitrary  magnitude-squared  characteristic  which  Is  approximated  optimally 
In  a  weighted  Chebyshev  (mlnimax)  sense.  Here,  the  optimization  procedure  can 
be  formulated  as  a  linear  programming  problem  and  the  filter  coefficients  can 
be  calculated  using  the  two-dimensional  discrete  Hilbert  transform  to  approxi¬ 
mately  factor  the  two-dimensional  magnitude-square  frequency  response  [4] . 

There  are  two  problems  in  using  this  algorithm.  The  first  Is  that  fairly 
large  amounts  of  computer  time  are  needed  to  design  filters.  Secondly,  the 
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factorization  of  magnitude-square  function  to  obtain  a  filter  Implementation  con¬ 
tains  an  Infinite  number  of  terms  and  hence  truncation  Is  necessary. 

Generally,  the  nonlinear  optimization  techniques  require  that  the  filter 
consist  of  first  and  second  order  sections  connected  in  cascade.  This  cascade 
form  has  many  advantages  over  the  direct  form  [1],  In  our  research  we  are  in¬ 
vestigating  the  approximation  problem,  i.e.,  how  well  can  a  cascade  of  first 
and  second  order  sections  approximate  an  arbitrary  frequency  response. 

Another  design  approach  is  to  find  filter  coefficients  which  approximate 
a  desired  frequency  response  without  adding  the  stability  constraint;  then,  if 
the  filter  is  unstable,  stabilize  it  so  that  the  resulting  filter  has  approxi¬ 
mately  the  same  frequency  response.  Several  stabilization  methods  exist.  The 
first  was  proposed  by  Shanks  [5].  It  stabilizes  an  unstable  filter  by  finding 
Its  double  planar  least  squares  inverse  (PLSl).  Until  recently  the  PLSI  of  a 
filter  was  assumed  to  be  always  stable.  However,  a  special  counter  example  has 
been  constructec'  by  Genin  and  Kamp  where  they  find  a  2x2  PLSI  of  a  4x4  array 
[6],  It  is  still  an  open  question  whether  a  PLSI  of  the  same  size  as  the  input 
array  Is  stable.  Jury  [7]  has  proved  this  for  the  special  cases  of  3x2  and  2x3. 
The  double  PLSI  of  an  array  will  produce  a  filter  array  whose  magnitude  response 
approximates  that  of  the  unstable  filter.  Howe/er,  in  many  applications  the 
approximation  is  not  adequate.  In  our  work  we  have  found  an  example  where  the 
double  PLSI  calculated  by  Read  and  Treitel  is  unstable.  But  no  claims  on  find¬ 
ing  a  significant  counter  example  will  be  made  until  Read  and  Treitel's  [8] 
calculation  of  PLSI  can  be  checked. 

Another  approach  stabilizes  an  array  by  developing  a  two-dimensional  dis¬ 
crete  Hilbert  transform  to  calculate  the  analytic  phase  function  from  the  mag¬ 
nitude-squared  frequency  response.  The  method  accomplishes  stabilization  with 
little  accompanying  distortion  of  its  amplitude  spectrum.  Several  developments 
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of  the  discrete  Hilbert  transform  exist.  Read  and  Treltel  [8]  generalized  It 
to  two  dimensions  by  a  direct  extension  of  the  one  dimensional  discrete  case. 
The  amplitude  spectrum  of  the  stabilized  filter  more  closely  approximates, the 
desired  magnitude  than  the  double  PLSI  filter.  However,  not  all  filters  pro¬ 
duced  by  this  method  are  stable.  An  example  of  this  was  given  in  their  paper 
and  we  have  found  many  more.  Another  method  relies  on  discretizing  the  con¬ 
tinuous  two-dimensional  Hilbert  transform  to  obtain  a  two-dimensional  discrete 
Hilbert  transform.  This  has  been  derived  by  Dudgeon  [9]  and  we  have  formulated 
and  programmed  a  special  case.  Preliminary  results  seem  to  indicate  that  this 
method  does  produce  stable  filters,  but  the  approximated  amplitude  is  notable 
distorted. 

We  have  been  studying  a  third  method  which  either  employs  the  cepstrum  or 
complex  cepstrum.  The  study  was  motivated  by  work  done  by  Ekstrom  and  Woods 
[10]  on  two  dimensional  spectral  factorization.  The  Inverse  cepstrum  of  the 
first  quardrant  of  the  autocorrelated  unstable  filter's  cepstrum  is  windowed 
to  give  a  stable  filter  whose  frequency  response  is  close  to  that  of  the  un¬ 
stable  filter.  Preliminary  results  show  that  even  though  stability  is  not 
guaranteed  the  resulting  filters  are  usually  stable.  This  method  can  also  be 
used  to  test  stability  of  any  two-dimensional  recursive  filter.  The  original 
filter  Is  stable  if  the  calculated  array  is  equal  to  the  original  array.  Be¬ 
cause  the  FFT  is  used  the  correspondence  is  only  approximate  and  some  equality 
measurement  must  be  used  to  ascertain  stability.  Through  experimentation  we 
found  that  for  3x3  arrays  the  filter  is  stable  If  the  mean  square  difference 
between  Input  and  output  arrays  is  less  than  .00000*f3»  (16x16  FFT  were  used.) 

Another  possible  way  of  checking  stability  and  stabilizing  unstable 
filters  is  to  work  with  the  complex  [12]  centrum  of  the  filter  array.  By 
properly  processing  the  cepstrum  we  can  obtain  arrays  which  are  nonzero 
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In  certain  distinguished  regions  In  the  (Zj„  Z^)  plane,  and  thus,  depending  on  the 
the  type  of  filter,  guaranteed  to  be  stability.  Many  problems  exist  In  one- 
dtmenslonal  cepstral  analysis  and  they  are  Increased  in  two  dimensions.  A  two 
dimensional  complex  cepstrum  program  has  been  written  using  a  modification  of 
a  new  phase  unwrapping  algorithm  which  was  developed  by  Tribolet  [11].  Details 
about  the  theory  and  Implementation  of  two-dimensional  cepstrial  analysis  will 
be  given  in  the  next  report. 

Future  work  will  include  a  detailed  analysis  of  two-dimensional  cepstral 
techniques  applied  to  filter  design  and  stabilization. 
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COMPARISON  OF  THE  PROJECTION  METHOD  WITH  SINGULAR 
VALUE  DECOMPOSITION 

S.P.  Berger  and  T.S.  Huang 

The  purpose  of  the  work  has  been  to  evaluate  the  relative  merits  of  the 
projection  method  of  image  restoration  as  compared  with  the  singular  value  de¬ 
composition  (SVD)  approach. 

The  projection  algorithm  is  an  Iterative  method  of  solving  a  set  of  linear 
equations,  where  the  equations  are  the  discrete  representation  of  the  degrada¬ 
tion  process.  The  algorithm  can  incorporate  a  priori  Information.  It  has  been 
shown  to  yield  effective  results  for  various  types  of  degradation.  Including 
space-variant  distortions. 

The  SVD  approach  has  been  used  successfully  in  image  restoration.  It  in¬ 
volves  the  treatment  of  the  degradation  In  matrix  form.  By  application  of  a 
pseudo- Inverse  matrix,  the  restoration  Is  hopefully  achieved. 

Both  approaches  are  limited  by  the  effects  of  noise.  The  required  number 
of  iterations  in  the  projection  method  and  the  number  of  terms  utilized  with 
SVD,  are  determined  subjectively.  The  effects  of  noise  Increase  with  the 
number  of  terms  and  the  number  of  iterations,  and  tend  to  overshadow  the 
restoration  process. 

In  the  actual  Implementation  of  the  SVD,  difficulties  with  this  approach 
have  arisen.  The  amount  of  computer  time  required  for  the  calculation  of 
eigenvalues  and  eigenvectors  can  be  prohibitive  for  large  degradation  matrices. 
Also,  perhaps  to  the  size  of  the  matrix,  the  actual  implementation  has  yielded 
faulty  results.  The  success  of  the  method  depends  heavily  on  the  type  of 
degradation  that  Is  effected  on  the  border  of  the  Image.  Minor  changes  in  the 
form  of  the  degradation  matrix  seem  to  create  major  problems  in  the  operation 
of  the  computer  implementation.  The  next  report  will  contain  at  least  a 
partial  resolution  of  these  difficulties. 
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FOURIER  DESCRIPTORS 
T.  Wallace  and  P.A.  WIntz 

The  Fourier  descriptor  (FD)  Is  one  method  of  describing  the  shape  of  a 
planar  figure.  Given  a  figure  In  the  complex  plane,  the  contour  can  be  traced, 
yielding  a  complex  function  of  time.  If  the  contour  Is  traced  repeatedly,  the 
periodic  function  which  results  can  be  expressed  In  a  Fourier  series.  Granlund 
[1]  defines  the  FD  of  a  contour  as  the  coefficients  of  this  Fourier  series. 

To  Implement  this  method  of  shape  description.  It  Is  necessary  to  sample 
the  contour  at  a  finite  number  of  points.  Since  the  discrete  Fourier  transform 
of  a  sequence  gives  us  the  values  of  the  Fourier  series  coefficients  of  the 
sequence,  assuming  It  to  be  periodic,  using  an  FFT  algorithm  satisfies  the 
definition  above.  The  computational  advantages  of  the  FFT  are  well  known. 

The  goal  of  this  work  Is  to  classify  the  shapes  of  objects  using  their 
Fourier  descriptors.  The  operations  of  rotation,  scaling,  and  moving  the 
starting  point  are  easily  implemented  In  the  frequency  domain  by  simple  arith¬ 
metic  on  the  frequency  domain  coefficients.  While  shapes  may  be  compared  In 
the  space  domain,  the  procedures  required  to  adjust  their  size  and  orientation 
are  computationally  very  expensive.  Normally  an  Iterative  type  of  algorithm  Is 
employed,  which  searches  for  an  optimum  match  Between  the  unknown  shape  and  the 
test  set. 

The  goal  of  classification  using  Fourier  descriptors  is  to  develop  an 
algorithm  which  will  normalize  the  size  and  orientation  of  a  shape  before  any 
comparisons  to  test  shapes  are  made.  If  this  can  be  accomplished,  the  classi¬ 
fication  process  becomes  a  simple  clustering  problem  with  no  Iterative  searches 
to  contend  with. 

In  our  last  quarterly  report,  the  FD  normalization  problem  was  dis¬ 
cussed,  and  the  results  of  an  experiment  presented.  It  was  shown  that  If  the 


contours  under  Investigation  were  bilaterally  symmetric,  much  less  information 
was  contained  In  the  phases  of  the  FD  than  In  the  magnitudes.  The  validity  of 
classifying  contours  using  the  magnitudes  of  their  FD’s  was  demonstrated  by 
classifying  18  airplane  contours  using  this  method.  The  Euclidean  distance 
measure  was  used,  and  the  distances  between  FDs  proved  to  correlate  well  with 
the  actual  differences  between  the  contours. 

The  goal  of  this  study  of  Fourier  descriptors  is  to  eventually  apply  this 
method  to  contours  traced  using  actual  photographic  data.  To  simulate  this 
situation,  an  experiment  was  conducted  using  20  aircraft  contours  digitized  to 
two  different  resolutions,  128x128,  and  64x64.  The  high  resolution  versions 
were  taken  to  be  accurate  representations,  and  the  lower  resolution  versions 
were  assumed  to  be  corrupted  by  noise  and  quantization  error.  Examination  of 
representative  contours  (Figs.  1-6)  show  that  the  128x128  contours  are  quite 
good  representations,  while  the  64x64  show  significant  distortion  of  the 
smaller  important  features.  This  experiment  was  performed  in  order  to  test 
various  distance  measures,  as  well  as  to  test  suitability  of  the  algorithm  for 
use  with  actual  photographic  data. 

While  the  experiment  described  in  [1]  was  useful  in  demonstrating  the 
general  validity  of  classifying  contours  using  distances  between  normalized 
Fourier  descriptors  (NFDS) ,  a  comparison  of  various  distance  measures  was 
difficult  due  to  lack  of  a  definitive  measure  of  similarity  among  the  contours 
themselves.  This  obstacle  was  overcome  by  using  the  different  resolution 
versions  of  the  same  planes. 

The  mean  square  criterion  used  previously  was  compared  to  an  absolute 
value  criterion.  The  results  using  the  absolute  value  criterion  were  slightly 
better,  as  every  64x64  contour  was  correctly  identified,  whereas  using  the 
mean  square  criterion,  only  19  out  of  20  were  correctly  identified. 


A  simple  application  of  Parseval's  theorem  shows  that  the  mean  square 
criterion  In  the  frequency  domain  corresponds  to  a  sum  of  the  point  by  point 
mean  square  error  in  the  time  domain.  There  are  some  difficulties  In  mathe¬ 
matically  relating  the  frequency  domain  absolute  value  criterion  to  the  con¬ 
tours  In  the  time  domain.  However,  it  Is  easy  to  compare  the  two  criteria  on 
a  qualitative  basis. 

It  Is  obvious  that  the  absolute  value  measure  will  be  more  tolerant  of 
one  or  two  large  coefficient  differences  than  the  mean  square  measure.  Using 
either  method,  the  largest  coefficients  will  account  for  most  of  the  distance 
measured.  Accordingly,  we  should  consider  the  possibility  of  large  variations 
between  coefficients  of  large  expected  value,  which  turn  out  to  be  the  lowest 
frequency  ones. 

The  second  largest  coefficient  for  each  airplane  FD  Is  A (—3) .  The  effects 
of  varying  this  coefficient  are  to  change  the  width  of  both  the  wings  and  the 
body  of  the  plane.  Smaller  detail  such  as  engine  shape  is  virtually  unaffect¬ 
ed.  Another  coefficient  of  large  expected  value  is  A(-1).  As  discussed  In 
[1],  this  coefficient  describes  the  length  to  wingspan  ratio  of  an  airplane 
contour.  More  generally,  varying  A (— 1 )  tends  to  elongate  any  contour.  Again, 
the  smaller  detail  Is  not  greatly  changed. 

In  summary,  the  absolute  value  criterion  should  tolerate  more  differences 
In  gross  structure  of  the  contour,  such  as  elongatedness  or  thickness  to 
length  ratio,  while  emphasizing  variations  in  smaller  detail.  It  seems  likely 
that  the  absolute  value  measure  might  correspond  more  clearly  to  the  differ¬ 
ences  which  human  observers  find  Important  than  does  the  mean  square  criterion. 

The  classifications  using  the  two  methods  were  too  similar  to  offer  a 
definitive  comparison  of  the  two  classification  methods.  However,  the  success 
of  the  experiment  shows  that  we  are  ready  to  apply  the  FD  algorithm  to  contours 
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extracted  from  actual  photographic  data.  We  plan  to  Interface  the  BLOB  con¬ 
tour  tracing  algorithm  to  the  FD  program  as  the  next  step  In  this  research. 
The  magnitude  vs.  phase  Information  question  will  also  be  examined  In  the 
context  of  developing  a  normalization  procedure  which  preserves  all  of  the 
Information  contained  In  the  contour. 
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FILTERING  TO  REMOVE  CLOUD  COVER  IN  SATELLITE  IMAGERY 
O*  R.  Mitchell  and  P,  L.  Chen 

I.  INTRODUCTION 

This  recent  work  Is  a  continuation  of  the  previous  report  "Satellite 
Imagery  Noise  Remov&l"  [1].  We  are  using  3-d imensional  homomorphic  filtering 
techniques  to  remove  cloud  cover  In  LANDSAT  data. 

In  the  previous  filtering  results  the  noise  power  spectrum  was  estimated 
by  classifying  each  region  of  the  noisy  picture  according  to  the  level  of  noise 
present  using  the  multispectral  data  analysis  software  system  developed  by 
Lars.  These  filtering  results  were  not  satisfactory  due  to  the  spectral 
indistinctabi 1 ] ty  of  clouds  and  concretes. 

It  may  also  prove  Impractical  to  use  a  multispectral  classification  pro¬ 
gram  to  find  the  noise  statistics.  Instead  it  may  be  possible  to  use  a 
generalized  cloud  power  spectrum  derived  by  averaging  many  sample  spectrums 
together. 

II.  INDIRECT  ESTIMATION  OF  NOISE  STATISTICS 

The  easiest  way  to  estimate  the  general  power  spectrum  of  cloud  is  done 
by  using  data  taken  over  water  where  the  reflection  is  almost  constant.  In 
this  case,  thi  transformed  scanner  image  (L  is  sun  illumination,  t  is  cloud 
transmission,  r  Is  ground  reflection,  and  s  is  the  received  scanner  Image) 

Log  [L  -  s(x,y) ]  «  Log  L  +  Log  t(x,y)  +  Log  [L  -  ar(x,y)] 

fs  reduced  to 

!  >g  [L  -  s(x,y)]  -  Log  [L  -  ar(x,y)]  +  K 

where  K  is  a  constant,  and  the  power  spectrum  obtained  is  that  of  noise  except 
for  the  d.c.  (0,0)  frequency  point. 
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This  power  spectrum  should  be  circularly  symmetric  since  clouds  have  no 
preferred  orientation  and  should  consist  mainly  of  low  spatial  frequency  com¬ 
ponents  since  clouds  are  relatively  large  and  smooth  functions  compared  to 
ground  reflectance. 

III.  THREE  DIMENSIONAL  FILTERING 

The  real  potential  In  the  cloud  filtering  process  is  incorporating  a 
third  dimension,  the  spectral  channels,  forming  a  three  dimensional  reflection 
r(x,y,z)  and  cloud  transmission  t(x,y,z).  The  generalized  linear  filter  thus 
employed  Is  three  dimensional  H(p,v,p)  using  three  frequencies  (two  spatial 
and  one  spectral).  Although  there  are  only  four  points  in  the  spectral 
dimension  for  LANDSAT  data,  the  method  has  good  promise,  because  most  clouds 
follow  a  fixed  response  in  the  spectral  dimension:  cloud  transmission 
Increases  with  wavelength  in  a  predictable  fashion.  When  this  information  is 
Incorporated  into  the  filter  (by  means  of  the  3-D  power  spectrum)  Image 
variations  which  have  the  cloud  spectral  response  are  filtered  out  and  image 
variations  which  do  not  follow  the  expected  response  of  clouds  in  the  spectral 
dimension  are  left  In.  The  three  dimensional  filter,  therefore,  tends  to 
reject  all  variations  that  are  low  frequency  in  the  spatial  dimension  and 
follow  the  cloud  spectral  response  In  the  third  dimension. 

Figures  1  to  I*  are  the  3-D  power  spectrum  obtained  from  averaging  the 
spectrum  from  three  separate  6Ax64  regions  of  clouds  over  water.  The  ordinate 
Is  a  log  scale.  Figure  I  represents  the  spectral  d.c.  slice  (sum  of  all 
four  spectral  points),  Fig.  2  represents  the  "one  cycle  per  spectral  band" 
slice,  etc.  One  interesting  observation  we  have  made  is  that  most  of  the 
cloud  information  Is  contained  in  slice  1  only  and  that  information  in  slices 
2,  3»  and  4  (other  than  the  d.c.  point  In  each)  represent  ground  reflectance 
Information  and  should  be  left  In  the  filtered  output. 
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Figures  5  to  8  are  the  3”D  filter  functions  obtained  for  one  particular 
64x64  section  of  the  Image  taken  over  land.  The  general  cloud  spectrum  was 
normalized  so  that  It  fell  below  the  spectrum  of  "signal  plus  noise"  and 
was  then  subtracted  from  It  for  the  signal  spectrum  estimate,  Iwte  the  filter 
tends  to  attenuate  low  frequencies  in  slice  one  but  leaves  them  In  slices  2 
through  4, 

Figure  9  shows  the  composite  filtered  for  21  64x64  blocks  of  LANDSAT  data 
scanned  over  the  Chicago  area  on  a  somewhat  cloudy  day. 

IV,  CONCLUSIONS 

Three  dimensional  filtering  of  multlspectral  data  to  remove  light  cloud 
cover  Is  a  distinct  possibility. 

In  Fig.  9  the  filtered  picture  shows  more  detail  than  the  original  cloudy 
one.  Lake  boundaries  and  highways  have  clearer  appearances  In  the  filtered 
pictures. 

The  model  of  the  cloud  distortion  needs  to  be  refined  based  on  the 
results  of  filtering  using  the  simple  model  presented  In  the  previous  report. 
It  may  be  necessary  to  consider  convolutional  effects  of  cloud  cover  as  well 
as  multiplicative  effects.  The  change  In  multlspectral  classification 
accuracy  after  filtering  may  be  a  suitable  measure  of  the  performance  of  such 
homomorphic  filter. 


Figure  I  3-0  cloud  power  spectrum  -  slice 
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Figure  3  *D  cloud  power  spectrum  -  slice 


Figure  4  3"D  cloud  power  spectrum  -  slice 
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Figure  9(a)  (Left)  3-D  filtered  output  (channel  3) 
(Right)  Original  Landsat  data  (channel 


Figure  9(b)  (Left)  3-D  filtered  output  (channel  4) 
(Right)  Original  Landsat  data  (channel 
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QTY  Manufacturer 

3  Beehive  Elect. 

2  Tax.  Inst. 

1  Dtgl-Data 


DEC 

DEC 

Fabrl tek 

Versa tek 
Comtal 

Data  Printer 
True-Data 
Tektronix 
DEC 


FACILITIES 

Description 

"Super-Bee"  Terminals 

"Silent  700"  Terminals 

Industry  standard  magnetic  tape  system; 

2,  9“track  and  1,  7“track  drives;  one  each 
NR2I  and  phase-encoded  formatters/controllers 

Dual -drive  DEC tape  unit 

RP03  disk  drive  (40  million  characters) 

96K-word  auxiliary  memory  system 
(64K  bought  by  ARPA,  32K  by  NASA) 

Electrostatic  matrix  printer 

Color  picture  display 

132  column,  600  L„P.M.  line  printer 

Punched  card  reader 

Model  4010,  graphics  display 

PDP-ll/45  computer  system;  system  includes: 
32K  memory 

FPP-ll  floating  point  processor  (NSF  money) 

H960  extension  mounting  cabinet 

3  ■  small  peripheral  mountings  blocks  (DD-ll) 

1  UNI  BUS  repeater/expander 

DHll,  16-1 ine  terminal  multiplexor 

KWll-p  programmed  clock 

"ANTS"  -  type  PDP-11/ IMP  interface 


Note:  Our  PDP-ll/45  Is  currently  operating  under  the  UNIX  system 
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